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ABSTRACT 



Context. We present a scaling argument that we develop into a shallow water theory of non-axisymmetric disturbances 
in annular sections of thin Keplerian disks. 

Aims. We develop a theoretical approach to understand physically the relationship between two-dimensional vortex 
dynamics that is known and their three-dimensional counterparts in Keplerian disks. 

Methods. Using asymptotic scaling arguments, varicose disturbances of a Keplerian disk are considered on radial and 
vertical scales consistent with the height of the disk while the azimuthal scales are the full 2n angular extent of the disk. 
For simplicity perturbations are assumed to be homentropic according to a polytropic equation of state. The timescales 
considered are long compared to the local disk rotation time. 

Results. The scalings relate to dynamics that are radially geostrophic and vertically hydrostatic. A potential vorticity 
quantity emerges and is shown to be conserved in a Lagrangian sense. Uniform potential vorticity linear solutions 
are explored and the theory is shown to contain an incarnation of the strato-rotational instability under channel 
flow conditions. Linearized solutions of a single defect on an infinite domain are developed and shown to support a 
propagating Rossby edgewave. Linear non-uniform potential vorticity solutions are also developed and shown to be 
similar in some respects to the dynamics of strictly two-dimensional inviscid flows. The relationship of the scalings and 
some of the resulting dynamics are considered with respect to other approximations employed in the literature. Based 
on the framework of this theory, arguments based on geophysical notions are presented to support the assertion that 
the strato-rotational instability is in a generic class of barotropic/baroclinic potential vorticity instabilities. Extensions 
of this formalism are also proposed. 

Conclusions. The shallow water formulation achieved by the asymptotic theory developed here opens a new approach 
in studying disk dynamics. 

Key words, accretion, accretion disks - instabilities 



1. Introduction 

The question of the formation of coherent structures in 
Keplerian disks has received increased attention in the past 
few years. Their importance has been recognized in playing 
a role in a number of processes including angular momen- 
tum transport and the accumulation of dust in the process 
of planet building. These coherent structures are persistent 
vortices and understanding their generation, structure, dy- 
namics and lifetimes has been the focus of a number of 
studies both in 2D (e.g. Godon & Livio, 1999, Bracco et al., 
1999, Umurhan & Regev, 2004, Johnson & Gammie, 2005, 
Petersen et al. A, 2007, Petersen et al. B, 2007, Bodo et al. 
2007, Lithwick, 2007 A) and recent 3D studies (Barranco 
& Marcus, 2005, Lithwick, 2007 B). 

There are many facets to the nature of vortices in 
strongly sheared and rotating environments like astrophysi- 
cal disks. One thing that most two-dimensional studies have 
shown, and which are in qualitative agreement with each 
other, is that from a patch of noisy initial conditions anti- 
cyclonic vortices naturally emerge and these structures are 
usually the sole occupants of radially demarcated zones. 
This generic phenomenon is observed in 2D shearing box 
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simulations (Umurhan & Regev, 2004, Johnson & Gammie, 
2005 Lithwick, 2007 A) and 2D cylindrical studies as well 
(Godon & Livio, 1999, Bracco et al., 1999, Petersen et al. 
A-B, 2007, Bodo et al. 2007). In all of these examples a vor- 
ticity (or potential vorticity) quantity primarily character- 
izes the flow dynamics. This vorticity may evolve (nearly) 
conservatively in a Lagrangian sense (e.g. Godon & Livio, 
1999, Bracco et al. 2000, Umurhan & Regev, 2004, Lithwick, 
2007 A) or may have a source term due to baroclinicity lead- 
ing to, for example, Rossby wave-like dynamics (i.e. wave 
generation due to radial gradients in temperature or en- 
tropy as in Lovelace et al., 1999, Li et al. 2000, Klahr & 
Bodenheimer, 2003, Petersen et al. A-B, 2007), or other 
wave generation processes of some type (see e.g. the com- 
pressible studies of Johnson & Gammie 2005, Bodo et al., 
2007). 

The dynamics of vortices in 3D simulations is less clear 
and the picture of what is occurring under these circum- 
stances is only now emerging. Lithwick (2007 B) has demon- 
strated that columnar anticyclonic vortices can manifest 
themselves in 3D shearing box simulations with no exter- 
nal gravity. It was also shown in that study that the sur- 
vival of the vortex column depends on the height of the 
column itself - if it is too large secondary instabilities may 
develop causing a toppling over of the structure, otherwise, 
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if it is short enough the structure continues to persist. In 
an anelastic stratified fluid model evolved in a shearing box 
with the presence of an external disk-vertical component 
of gravity, Barranco & Marcus (2005) show that columnar 
anticyclonic vortices show a tendency to breakdown but are 
replaced with anticyclonic structures of restricted vertical 
extent located away from the disk midplane. The way that 
this process emerges is subtle and a full mechanical under- 
standing of it still awaits - perhaps aided by a shallow water 
model description. 

Shallow water models are effective at capturing salient 
physical effects involving stratification without having to 
resort to full 3D analyses and simulations. In this way 
they have proven to be indispensiblc tools in furthering our 
physical understanding of oceanic and atmospheric flows 
(Pedlosky, 1987) and have recently been used in studies of 
cxoplanet atmospheres (Langton & Laughlin, 2007, Cho et 
al., 2007). 

In a series of studies (Balmforth ct al. 1992, Balmforth 
et al., 1993, Balmforth et al. ,1995, Balmforth & Spiegel, 
1996) a shallow water formalism was presented and ana- 
lyzed to study the axisymmetric structure and dynamics of 
thin cold disks both with and without the effects of self- 
gravity. This formalism proved to be especially effective at 
capturing and describing the dynamical response of self- 
gravitating disk systems. In a limit akin to a shallow wa- 
ter reduction, the linearized equations of motion of small 
sections of accretion tori have been shown to lead to inter- 
esting non-axisymmetric dynamics including processes such 
as the Papaloizou-Pringle instability (PPI, Papalouizou & 
Pringle, 1984,1985, Goldreich et al., 1986). We shall return 
to this in the Discussion section. 

Owing to the ubiquity of (potential) vortex dynamics 
in 2D and the appearance of long-lived coherent structures 
in 3D simulations it is worthwhile to develop a framework 
to better understand the formation of vortices in 3D sce- 
narios by comparison to the processes leading to 2D vor- 
tex structures. A shallow water theory appropriate for non- 
axisymmetric disturbances of an accretion disk should serve 
a critical intermediary role in understanding the physical 
development of such processes. In this context, meteorolog- 
ical studies ought not to be overlooked: potential vorticity 
dynamics ( sometimes referred to in this work as "PV") 
and its role in the generation and development of plane- 
tary storms in the upper atmosphere have been long ap- 
preciated (Hoskins et al., 1985, Bishop & Davies, 1994, and 
many others). 

The road to such a shallow water formalism is guided 
in part by the results thus far established in the literature. 
Especially important is the tendency seen in 2D cylindri- 
cal studies for (anticyclonic) vortices to either fill out or be 
the sole occupant of the azimuthal scale of their radially 
restricted domain (unambiguous examples of this may be 
found in Godon & Livio, 1999, Petersen et al. A, 2007, and 
Petersen et al. B, 2007). When viewed in a reference frame 
moving with a steady vortex the scale of the azimuthal 
velocities ("£/") are larger than the corresponding radial 
velocities in proportion to the ratio of the azimuthal to ra- 
dial length scales of the vortex structure. Assuming that 
the azimuthal scale is the radial position where the vortex 
is found ("i?o") then the circulation time of these vortex 
structures is ~ Ro/U. Such vortices in thin (cold) disks, 
whose azimuthal velocity scales are no greater than the typ- 



ical sound speeds, will have time scales which are 1/e longer 
than the typical rotation time of the disk at the location 
of the vortex, where e is the ratio of the local soundspeed 
to the local Keplerian velocity. In typical cold disks this 
number is on the order of 1/20. 

We therefore reanalyze the full equations of motion (in 
cylindrical coordinates) by considering the dynamics tak- 
ing place in an annular section of the disk instead of the 
usual box section known as the Shearing Box (Goldreich & 
Lynden-Bell,1965). The radial and vertical scales are scaled 
by the disk height, H n = ei? 0l where (as above) Rq is the 
radial scale of the disk. The azimuthal length scale consid- 
ered is the full 2ttRo complement of the disk. The veloc- 
ity perturbations sitting atop the basic flow (which we as- 
sume to be Keplerian, i.e. rotationally supported) are such 
that the azimuthal velocity fluctuations scale with the lo- 
cal soundspeed while the radial and vertical velocities are 
much smaller than this by order e. Most importantly, tem- 
poral variations do not scale by the usual rotation time 
scale (given by the inverse of the Keplerian rotation fre- 
quency flo °f the disk at R ) but on a time scale which 
is 1/e longer. In order to affect a shallow water reduction 
in a transparent way without entirely dispensing with the 
compressible nature of the gas itself we assume here that 
the disk gas behaves homentropically where, in particular, 
we assume a polytropic equation of state characterized by 
the index 7 (as in Salby 1989, Cho et al. 2007 and see also 
the earlier studies of Papalouizou & Pringle, 1984, 1985, 
Goldreich et al. 1986). 

These adopted scalings, especially the slow timescale, 
lead to conditions in which disturbances are vertically hy- 
drostatic. It also follows that disturbances are geostrophic 
but only in the radial direction. This means to say that the 
dynamics are everywhere in radial geostrophic balance, i.e. 
radial pressure gradients balancing Coriolis effects. These 
circumstances are similar to the situation encountered re- 
sulting from similar scalings of the standard equations of 
the shearing box (e.g. the Large-Shearing Box limit in 
Umurhan & Regev, 2004, Umurhan, 2006). This type of 
leading order radial geostrophic balance has been encoun- 
tered in prior related work as well (Balmforth & Spiegel, 
1996, Barranco et al., 2000). 

The main result of this analysis then shows that in this 
asymptotic limit there exists a potential vorticity quantity 
which is conserved in a Lagrangian sense in the radial and 
azimuthal directions. In other words we have 

- _ n (2 - g) + d x v 
Dt U ' " h & W 

where 5 is the potential vorticity. The azimuthal velocity 
v represents perturbations atop the basic (local) Keplerian 
flow and where h(x, y, t) is the height of the disk measured 
from the disk midplane. The Lagrangian time derivative 
operation is restricted to the two horizontal directions: 

— = d t + ud x + {v- qQ, x)d y , 

and with fi measuring the local rotation rate of the disk 
(here it is scaled to 1). The basic Keplerian flow is given by 
— qQox in which x is the local (Cartesion) radial coordinate 
and g = 3/2. The azimuthal and radial velocities are given 
diagnostically in terms of the height of the disk, 

v = V(h) = ^d x h 2 , (2) 
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in which the operators V and T are 
T{h) = (d x V - qtl - 2^o) • dyQ h + 2d x 



(3) 



7-1 

7 + 



jQhdyV 



V = ^2(2 - q) - 2fl V ■ d x - 2d x 



7+1 



Qhd x ). (4) 



with Q h = \Q 2 h 2 . 

The radial and horizontal velocities are derived to be 
related to the vertical linearly varying velocity, (i.e. w — 
Q z z, where Vt z — £l z (x, y, t)) via a quasi-anelastic equation 
ol state, 



Or 



d y v + - ^ \ d z v = 0. 



7 



1 



(5) 



As such, these equations describe varicose disturbances 
with respect to the disk midplane (cf. usage in Balmforth 
et al., 1992) since the vertical velocity has odd symmetry 
with respect to the midplane. 

We use the set (JlJSJl to solve a number ol example 
problems. A linearized study ol uniform potential vorticity 
disturbances in a channel flow configuration is presented. 
Disturbances with only one solid boundary on a semi- 
infinite domain manifest as a single edgewave (Goldreich 
et al. 1986) propagating along the wall. Disturbances on 
a finite domain sandwiched by two solid boundaries man- 
ifest as two edgewaves - each associated with their re- 
spective boundary. These objects are similar to the edge- 
waves discussed by Goldreich et al. (1986) rationalizing 
the Papalouizou-Pringle instability (Papalouizou & Pringle, 
1985). These edgewaves become unstable if the width ol the 
channel lies between a narrow band of values less than the 
height of the fluid layer. For these circumstances we under- 
stand this instability as being the incarnation of the stra- 
torotational instability discussed by Dubrulle et al. (2004). 

Linearized disturbances of a radially infinite domain 
with a single global potential vorticity defect appear as a 
Rossby edgewave propagating along the defect. When the 
potential vorticity on either side of the defect becomes in- 
finite (negatively) the Rossby edgewave limits to the edge- 
wave solution associated with the semi-infinite domain dis- 
cussed above. We therefore rationalize the use of a solid 
boundary as representing this limiting form. 

We also investigate the dynamics of non-uniform poten- 
tial vorticity disturbances subject to the popular shearing- 
sheet boundary conditions (Goldreich & Lynden-Bcll, 
1965). We find that the potential vorticity evolves in ways 
that are very similar to 2D vorticity disturbances as in 
Chagelishvilli et al. (2003) and Umurhan & Regev (2004). 
We find that potential vorticity waves which are leading (in 
character) show strong transient growth thanks to the Orr- 
tilting mechanism (Schmid & Hcnningson, 2000) which is 
inline with the 2D studies cited as well as being inline with 
the general tendency for this to happen in 3D linearized dis- 
turbances as studied by Korycansky (1992) and Tevzadze 
et al. (2004). 

This work is organized as follows. In Section 2 we de- 
velop the derivation of (P(5]). In Sections 3-5 we present the 
linear results of the study of the uniform, single defect and 
non- uniform potential vorticity disturbances (respectively) . 
In Section 6 we, summarize our results, discuss the relation- 
ship of the scaling assumptions employed here in context 



of other approximations employed in the literature, con- 
sider some of the dynamical results with respect to other 
published results and we discuss possible future directions 
and extensions. We suggest that despite the use of artificial 
boundary conditions in the classical SRI analysis, that the 
instability could be possible in real disks supporting regions 
of oppositely signed potential vorticity gradients. 

2. Scaling Analysis and Derivation of the reduced 
equations 



The full equations of motion 
dp 



dT 



pV • U = 0, 

dU _ 1 
~dT ~ p 
dS 

dT=°> 



VP- V$, 



(6) 
(7) 
(8) 



are considered in cylindrical coordinates (R, </>, Z) where 
U = (Ur, U<p, Uz) and where the entropy S = Cy InP/p 7 
with 7 = Cp/Cy, i.e. the ratio of specific heats. The ex- 
ternal potential comes from the central gravitating source 
and has the form 



Q R 



R 2 , Z 2 



1/2 ' 



in which the frequency f2o = y GM/Rq is the Keplcrian 
rotation rate at the fiducial radius R = Rq. The derivation 
of the equations follows that of Umurhan & Regev (2004) . 
We shall consider a region centered around the radial point 
R = Ro and Z = 0. We transform our reference frame so 
that we are rotating with frequency Clo appropriate for the 
radial point Rq. We begin by stating a number of scalings 
for the various quantities involved. We assume that at this 
point the density of the fluid has a scale p and that the pres- 
sure scales as pc 2 s where c s is the sound speed. Furthermore, 
we define the non-dimensional parameter 



Vk ' 



v k = ii R , 



where v/- is the Keplerian rotation velocity at R = Rq. 
We shall assume that e is less than 1 and "small" . We shall 
consider radial and vertical extent which scale as ei?o - thus 
this is formally a box-like region in these directions. For the 
azimuthal scales we shall assume them to be of the order R . 
In this sense the domain represents a thin annulus instead 
of a classic small box like the shearing box. 

The departures from the classic shearing box equation 
derivation (e.g. as in Umurhan & Regev, 2004) appears at 
this point. Following the heuristic observations and argu- 
ments proffered in the Introduction, we assume that time 
scales as 1/eOo and that the vertical and radial velocities 
(in the frame moving with frequency S7o) scale as ec s while 
the azimuthal velocity fluctuations scale as the usual c s . In 
summary we assume the following scalings 

R — > Rq + eRqx, Z — » eRoz, 

0^y-n o t, T^t/(eQ ), 

U R — > ec s u, U z — > ec s w, — ► v k + c s v y , 



pp, P -> pc s p. 



(9) 
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The variables p, p,t,x, z,y,u,w,v y appearing on the 
righthandsides of all transformations are understood to be 
order 1 non-dimensional quantities corresponding to their 
dimensional counterparts. Inserting these expansions into 
the governing equations of motion expressed in the rotat- 
ing frame yields the following set of equations at lowest 
order in e, 



dp „ 

— 2Qq% 



+ O(e) 

—d x p + 2f&qx + O (e 2 ) 
p 



^ + 2Q u = -ld yP + 0(e 2 ) 
at p v ' 

= --d zP - fl 2 z + O (e 2 ) 
P 



(10) 
(11) 
(12) 

(13) 
(14) 



where s = lnp/p 7 , q = 3/2 and Oo = 1. The (Cartesian) 
time derivative operation is given to be 

— = d t + ud x + v y dy + wd z , 

while the divergence operation is 
V • u = d x u + dyVy + d z w. 

We shall assume that the fluid behaves as a barotrope 
where, in particular, we assume that it behaves like a poly- 
tropic gas, i.e. 



P = Kp\ 



(15) 



and consequently it means that the fluid layer is homen- 
tropic. This is to say that the fluid layer is characterized 
by a constant entropy given to be s = In if. We may now 
exactly solve for (fT3|) by defining the fluid enthalpy Q (as 
in Goldreich et al., 1986), 



dQ 



-dP = K-l—dpi- 1 . 

P 7-1 



In terms of the enthalpy Q13p is simply 
d z Q = -OqZ, 
yielding the solution 



i^^-z 2 ) 



(16) 



(17) 



(18) 



The height h = h(x,y,t) defines where the fluid density 
and pressure go to zero and is a function of the horizontal 
coordinates and time. 

The entropy equation (|14p . the equation of continuity 
(JTOj) , the polytropic gas assumption and the solution to Q 
(fl~8"|) , all taken together yield 



dQ 
dt 



( 7 -l)QV-u = 0. 



(19) 



By choosing this form for the vertical velocity we are re- 
stricting attention to varicose disturbances of this annular 
disk section. Explicitly inserting the solution form (|18[) into 
(|19|) and making use of the ansatz for w we find 



For the solutions we shall examine here, as done in the 
derivation of shallow water equations, we state the ansatz 
that the horizontal velocities have no z-dependencies while 
the vertical velocity has a linear z-dependence, i.e. 

u = u(x,y,t), v y = v y (x,y,t), w = Cl z (x,y, t)z. (20) 



D 
Dt 



( 7 -l)|fi 2 (/i 2 
in which 



Z 2 )(d X U + OyVy + Q x ) = 0, 



(21) 



D 
Dt 



= d t + ud x 



Vydy. 



On account of their linear independence, we equate to zero 
the coefficients of each power of z appearing in (|2"Tj) . In this 
particular situation this yields two equations, the first is 
the equation for the height h 



Dh 
151 



hVt z 



(22) 



while the second is an expression resembling an anelastic 
condition 



8 X U + OyVy + 



7 + 1 
7-1' 



n z = o. 



(23) 



We note that the limit 7 — ► 00 recovers the condition of 
incompressible fluid flow. Addressing (jll[) we rewrite the 
azimuthal velocity v y = —qQox + v to find this equation 
appearing now as the geostrophic balance discussed in the 
Introduction, 



- 2tt v = -d x Q, 
while (fl~2"|) now appears 
Dv 



Dt 



+ Oq(2 - q)u = -d y Q. 



(24) 



(25) 



Note that the time derivative looks like 

D 
Dt 



(dt — qfloxdy) + ud x + vd y 



The set of equations (|22M25j) form the basic equations of 
this asymptotic limit. 

We may proceed even further and develop a potential 
vorticity equation. We begin by operating on (|24|) with d y 
and subtracting from this the result of operating on (|25p 
with d x leaving 



Dd x 



Dt 



+ [fi (2 -q) + d x v] (d x u + d y v) = 0. 



(26) 



Using (|23p in the above and combining the result with 
finally results in fl} . The quantity 5 is the potential vortic- 
ity and |T]) is the analogous Taylor-Proudman theorem for 
this shallow water system (cf. Pedlosky, 1987, Cho et al., 
2007). We may further isolate the radial velocity in terms 
of h. We do this by (i) multiplying ([22]) by (ii) oper- 

ating on the resulting equation with the partial derivative 
operation d x , (iii) subtracting the resulting equation from 
([25]) , (iv) reordering the result and using (|24|) results in the 
diagnostic equations for u and v found in (PE]). 
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3. Constant uniform potential vorticity solutions 

We consider the dynamics emerging from a situation where 
5 is uniform and constant. We imagine that we disturb 
a uniform equilibrium h = ho- Thus this means that the 
uniform potential vorticity state is given by 



- -n 2 ' q 

h 7-1 



(27) 



Through the use of ([24]) we find that the equation to be 
solved is 



n (2 - 1) + n d 2 x (\h 2 



■y + l 

E h— = 0, 



(28) 



subject to specified boundary conditions. Given that we 
will be considering disturbances around the uniform state 
we define a new quantity O as 

h — hp, 



h 

Thus, the potential vorticity equation (|28|) is written more 
simply as 



hhl%(l + 0f+2{2-q) 1- (1 + 0)5 



0. 



(29) 



We shall consider here only two possible boundary con- 
ditions. The first is that the normal velocity is zero at the 
position x — and that the pressure remain fixed as x — > oo 
which translates to requiring h — > ho- We shall refer to this 
scenario as the "semi-infinite domain" . Referring to ([23]) the 
condition on the normal velocity is the same as requiring 



(d t - qn xd y )v + vd y v + \d v ttl(l + 9f = 0, 



(30) 



at x — and where v = ilojd x (l + 9) 2 . The second set of 
boundary condition is that we require the normal velocities 
be zero at the positions x = ±A. In other words, we require 
that ([30]) be satisfied at those boundaries. We refer to this 
as the "channel flow" . 

We consider linear disturbances for 9. The governing 
equation for the disturbance PV, S', becomes 



E' ee hldlO - 2(2 - q) 



7 + 1 
7-1 



= 0. 



(31) 



The normal velocity condition at the boundaries becomes 



(d t - qn xd y )d x 9 + 2n d y 9 = -2(2 - q)u = 0. 



(32) 



Semi-infinite domain. The solution appropriate to the 
semiinfinite domain is 



e 



2(2 -g) 7+1 

K 7-1' 



(33) 



where we select the k > solution. The amplitude time 
evolution of the amplitude A is obtained from the velocity 
boundary condition 



(34) 



d t A + 2n K- 1 d y A = 0. 



Making the normal mode ansatz of the form A(y, t) ~ 
expia(y — ct) where a is the streamwise wavenumber, we 
find that the wavespeed is given by 



-2VLo/k. 



(35) 



Since c < 0, it follows that these are wave patterns which 
propagate in the negative y direction, that is to say, against 
the mean flow. 

Channel flow. Under channel flow conditions the 
normal-mode solution of 9 is 



A cosh k(x + A) — B cosh k(x — A) e 



Enforcing the linear boundary conditions 
x = ±A using the solution ([36]) we find 



ia(y-ct) 



(36) 



at the points 



- (c + gO A) An sinh 2kA + 2f2 (-4 cosh 2kA - B) = 0, (37) 
at x — —A and 

- (c - qfl A)BK sinh 2nA + 2fl (A - B cosh 2kA) = 0, (38) 

at x = A. These two simultaneous equations are satisfied 
non-trivially if the wave speed c satisfies, 



q 2 QlA 2 H ^ ( 1 - qnA coth 2kA 



(39) 



The wavespeeds approach the local mean flow speed of 
the respective boundaries as the separation between the 
walls becomes large, i.e. 



lim c 2 = fllA 2 q 2 



(40) 



On the other hand, when the wall separation becomes small 
the wavespeeds limit to 



1 



lim c 2 = hnQ?, 

A^O °7 + l 



(41) 



Examining the wavespeed prescription (|39p in the case 
where there is no shear, i.e. q = 0, shows that there is 
no possibility of instability. However, we do see that there 
seems to exist a band of instability in the gap width A for 
values of q > 0. In FigurcQ]we plot the region of instability 
and growth rate as a function of a variety of parameters. 
We see that range in A for which instability is predicted is 
less than ho', for instance, for 7 = 5/3 instability is expected 
approximately only when 0.26/io < A < 0.37/io- The range 
widens as 7 increases but is still less than ho for physi- 
cally reasonable values of 7. Note that the thickness of the 
channel is measured by 2 A. 



4. Normal-modes of an infinite domain with a 
single PV defect 

In this section we consider the response of an infinite do- 
main separated by a single defect in the potential vorticity. 
For simplicity we again treat ho as constant. Therefore, the 
defect appears in this case as a jump in the basic state 
shear, i.e. 



q-, 



X > 0, 
x < 0, 



Vo 



-Qoq+x, x > 0, 



-a 



q~x, 



x < 0, 



(42) 



Where Vq is the steady state uniform shear. This formu- 
lation ensures that the background velocity is everywhere 
continuous. Normal mode disturbances can be developed 
for this configuration much in the manner developed in 
the previous section. This means separately solving (f3"Tj) 
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Contours of Growth Rates: lm[c/(h o n o )] 




0.05 0.1 



Fig. 1. Contours of growth rates of the SRI as a function 
of 7 and A/ho- The growth rates are scaled by fio^o- The 
behavior for q = 3/2 is shown. The shaded region corre- 
sponds to the instability band. Horizontal line designating 
7 = 5/3 is also shown. 

on either side of the defect and then matching the solu- 
tions across the interface. Q Thus, for x > 0, disturbances 
which decay as x — > oo are 



= A + e ia{v ~ ct) e 



2(2- 9+ ) 7 + 1 

7-r 



h 2 



(43) 



while for x < 0, disturbances decaying as x — > — oo are 

2 _ 2(2 7 + 1 



h 2 



7" 



1 



while in both cases the positive root for k± is chosen. Since 
9 represents the perturbed height of the disk we require 
that it be continuous across the interface - this results in 
A- = A + = A. Because we are investigating normal mode 
disturbances, it also makes physical sense to match the nor- 
mal velocities across the interface. There are a number of 
ways to do this but the mathematically simplest is to con- 
sider the linearized version of 11121) . 

(45) 

Since the linearized expression of the geostrophic balance 
relationship is v' = \d x 9, in order for the radial velocity u' 
to be continuous across the defect we must have that 



(dt — qfloxd y )v + flo(2 — q)u' = —h d y 9. 





x^0 + 


V 









ia(—c — qQox)d x 8 + 2Qqio9 



.0+ 



0. (46) 



Given the solution for 9 on either side of a; = 0, in order 
for this to be satisfied is if the wavespeed c satisfies the 
relationship 



-29. 



q+ - q- 



.(2 -?_) + K_(2-q + \ 



(47) 



1 Another way to look at this would be to argue that if the 
perturbed height of the disk were to change discontinuously then 
it would correspond to a delta-function horizontal pressure gra- 
dient. 



The wavespeed represents the propagation of a Rossby 
edgewave in the sense presented in Hoskins et al. (1985). 
So long as q± < 2 the edgewave propagates in a negative 
sense with respect to the coordinate y if q+ > q- . Stated in 
another way, when this condition is met the Rossby wave 
travels against the mean flow. 

We note a significant mathematical limit of the 
wavespeed prescription (|47|) : when g_ — > — oo, while keep- 
ing <7 + finite, we find that 



— 2f2rj/ K H 



(48) 



By identifying q + with q it follows from the definition of k 
P3p that k + = k. Hence, we have recovered the wavespeed 
form for the semi-infinite domain, i.e. (|35|) in which the 
edgewave propagates against the local mean flow. This view 
allows us to interpret the presence of an inner wall bound- 
ary condition to be like a PV-defect whose PV on the one 
side of the defect (e.g. x < 0) greatly exceeds in magnitude 
the PV of the other side of the defect. 

This implication can also be understood by studying 
(I46[) : when, say, — g_ 3> 1 then it means that the normal 
velocity approaching the defect from + should be corre- 
spondingly small. In the infinite limit it is equivalent to 
requiring u'(x — > + ) = 0. 



5. Linear, non-uniform potential vorticity solutions 

We consider here linearized solutions to the potential vor- 
ticity 5' when it is not in general zero. This means solving 



(d t - qflody)-' = 



(49) 



given the definition of 3' to be found in pip , this may be 
rewritten instead as 



(44) (dt - q9 xdy) 



2(2 - q)l±l - hldl 



7-1 

The radial velocity is the solution to 



= 0. 



(50) 



2(2-9)- 



7 



1 



h Slo 



7+1 
-2(9 + 2) 



hldl 



7 + 1 



(51) 



We shall consider here wave solutions on the classical 
sheared periodic domain as introduced by Goldreich & 
Lynden-Bell (1965). To this end we define a new coordi- 
nate system in which the new azimuthal coordinate moves 
with the local shear 

X = x, Y ~ y + q n a xt, T = t, (52) 
so that all derivative operations now translate to 



d t -> &r, 



d x + qrioTd v , 



d„ -> dy 



(53) 



Then in the usual way we find the evolution equation for 
H' can now be written as 



dn 



h 2 (d x + qn Td Y ) 2 + 2(2 -q) 



7+1 



9 = 0. 



(54) 



7-1. 

In this frame we assume Fourier solutions of the form 

9 = 6 H (T)e ikX+iiY , (55) 
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where k and I are wavevectors. The solution to the wave 
amplitude 9ke{T) is given simply by 



h\k 2 



2(2 



>kt 



'k( 



h§(Af + gno£T) 2 + 2(2-g)2±i' 



where 9ki is the wave amplitude at time T = 0. The solution 
to the aziumuthal velocity may also be expressed as 



= -ij£loho(k + qQ, Q tT)6ki- 
Similary, representing the solution for u as 
Uke(J- )e , 
we have from (1511) 



(57) 



Uk£ 



2(2 + q ) + -^h 2 (k + q n eTf 
2(2 - q ) + ^h 2 (k + q n ery 



'kl- 



(58) 



Long time evolution of the solution to 9ki and Uki shows 
decay as 1/T 2 while vm decays as 1/T. If one restricts at- 
tention to solutions where T > 0, leading waves, i.e. those 
in which k£ < 0, will show a maximum in the amplitude of 
9ki when r max = k/(qD,Qi). The ratio of the amplitudes of 
these trailing waves to their initial amplitude is 



7 kt 



(0) 



= 1 



h\k 2 



7' 



1 



2(2-g) 7 +r 



showing that amplification is greater as k increases. 
Furthermore, at r max we see that Vke = 0. 

On the other hand, trailing waves, i.e. those in 
which kl > 0, show monotonic decay for all T > 
0. Examples of this sort of behavior is well-known 
for the evolution of two dimensional vorticity perturba- 
tions (e.g. Chagelishvilli et al. 2003, Umurhan & Regev, 
2004) and three-dimensional disturbances (Knobloch, 1984, 
Korycansky, 1992, Sternberg, 2005, Balbus & Hawley, 
2006). 



6. Discussion and Summary 

In 2D inviscid, incompressible vortex flow in a background 
linear shear profile, the vorticity is a quantity which is 
conserved in a Lagrangian sense. In terms of the non- 
dimensional quantities we have come to use in this study, 
this means to say that that 



[dt + ud x + (v — qVtox)d y ] u> = 0, 



(59) 



in which the vorticity lo may be expressed in terms of a 
stream-function tp according to u = (d 2 + dy)if). The radial 
and azimuthal velocities are then given by u — d y ip, v y = 
—d x ip. By construction the u and v y velocities satisfy the 
criterion of incompressibility. 

By comparison, the asymptotic system we have devel- 
oped here, representing the varicose disturbances of an an- 
nular section of a cold rotationally balanced accretion disk 
containing a homcntropic gas, has a similar form. In other 
words it means to say that the corresponding Lagrangian 
conserved quantity is the potential vorticity S, i.e. 



[d t + ud x + (v - q n x)d y ] 3 = 0. 



(60) 



The potential vorticity is related to the disk height h ac- 
cording to 



(56) Q: 



n (2 -q) + d x V(h) 



2 + 1 

hi- 1 



(61) 



where 7 is the polytropic index characterizing the home- 
ntropic layer under consideration. The radial and vertical 
velocities are given by 



v = V(h), 



(62) 



in which the operators appearing above are given in (2- 
4). The disk height h is, in this case, analogous to the 
streamfunction in classical 2D inviscid incompressible the- 
ory. Unlike the latter, the horizontal velocities reflect the 
dynamics of a quasi-anelastic relationship of the flow quan- 
tities in 3D, i.e. 



7 + 1 

8 X U + OyVy ~\ ~fl Z = 0, 



7' 



(63) 



where in general the vertical velocity is given as w = Q z z. 
The relationships (|6H63[) including the definition of the op- 
erators (2-4) are self-contained. 

Note that the quantity f2 z does not play a direct role in 
the evolution of the potential vorticity. Its meaning can be 
understood by noting that the expressions establishing the 
relationships between the horizontal velocities and the disk 
height h are consistent with the material derivative of the 
disk height, that is, its role may be understood in terms of 
the relationship, 



[dt + ud x + (v — qQ,Qx)d v ] h = w(z = h) = hfl z . 



6.1. Relationship to other approximations 



(64) 



In Table 1 we present a comparison of the various approx- 
imations used to model small sections of accretion disks 
(i.e. conditions in which e < 1). The four approximations 
presented are ones used in the literature and we describe 
them briefly below. In all of the following standard approx- 
imations the vertical extent of the disk H is assumed to be 
smaller than the radial scales of the disk by order e - that is, 
H ~ 0(eRo). Additionally, and within the moving frame 
assumed throughout these considerations, a Mach number 
is defined by, 8 = U/c s , where U is the largest velocity scale 
encountered in the frame. 

• SSB. The small shearing box equations are the classi- 

cal ones used to investigate the dynamics on very small 
scale sections of disks. Originally proposed by Goldreich 
& Lynden-Bell (1965) they have been the basis of nu- 
merous studies pertaining to the question of hydrody- 
namic disk turbulence (Lesur & Longaretti, 2005) and 
nature and maintanence of 3D disk vortices (Lithwick, 
2007 B). A scaling analysis (Umurhan & Regev, 2004) 
shows that these disturbances are assumed to have 
Mach numbers smaller than the aspect ratio of disks 
(i.e. S <C e). The corresponding length scales of distur- 
bances are smaller than the vertical extent of the disk 
as well but the time scales of all disturbances are equal 
to the local disk rotation times (i.e. ~ l/f2o)- 

• LSB. The so-called large shearing box equations are 

derived from a scaling analysis similar to those lead- 
ing to the SSB (Umurhan & Regev, 2004) except that 
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the velocities are scaled to have Mach numbers which 
are order 1. The corresponding length scales of distur- 
bances compare to the vertical scale of the disk. As with 
the SSB, the time scales of disturbances equal the lo- 
cal disk rotation times. The resulting equations are the 
compressible version of the equations of the SSB which 
also allow for the effects of gravity waves as there is now 
also an entropy equation which is evolved as well. 

• Elongated Vortex. Assuming a flow Mach number which 

is small, an analysis of a Keplerian disk in which one 
considers disturbances whose azimuthal and vertical 
scales compare to the local disk height while the ra- 
dial length scales are smaller by an order 6 lead to the 
equations describing elongated vortices (Barranco, et 
al., 2000). Under these circumstances, while the sub- 
sonic horizontal and vertical velocities of disturbances 
scale by the Mach number of the soundspeed (i.e. Sc s ) 
the scale of the radial velocity disturbances are smaller 
than these by another order of the Mach number. 
Consequently the time scale of all disturbances are an 
order 1/6 longer than the local disk rotation time. 

• Anelastic Model. Barranco & Marcus (2005) present a 

scaling argument leading to a set equations which we 
refer to here as the Anelastic Model. Flow quantities 
are formally assumed to have small Mach numbers (i.e. 
i5 < 1). Azimuthal and radial disturbance scales are as- 
sumed to be an order 5 smaller than the vertical scale 
of the disk H while the vertical disturbances scales are 
of the order H . The subsonic azimuthal and radial ve- 
locity disturbances scale by the Mach number factor of 
the soundspeed while the vertical velocities are a fac- 
tor of 5 even smaller than these. Consequently, the time 
scales of all these vortical disturbances is of the order of 
the disk rotation time. By assuming that the dynamical 
pressure and densities scale as they do one arrives at a 
circumstance in which sound waves are filtered from the 
equation set resulting in equations of motion describing 
an anelastic fluid. 

The scalings describing the dynamical response of the fluid 
pressure and density disturbances are also presented in the 
table for the purpose of comparison. The point to take away 
from this is that the pressure and density fluctuations are 
small by order S 2 their respective steady state values in the 
equations describing elongated vortices and the Anelastic 
Model as compared to the corresponding pressure and den- 
sity disturbances considered in this work and for the LSB. 
The dynamical pressure fluctuations also scale as a factor 
of 5 in the SSB. The small scalings of the thermodynamic 
quantities affect the filtering out of soundwaves. 

The scaling analysis we have employed here to arrive at 
this shallow water limit have some features similar to the 
arguments leading to the the Elongated Vortex equations of 
Barranco et al. (2000). Physically it means that there exists 
radial geostrophic balance of dynamical quantities - a con- 
dition which has been visited and explored in prior studies 
(Balmforth et al., 1992, Balmforth & Spiegel, 1996). The 
important feature that the approximations made in this 
study shares with the arguments leading to the Elongated 
Vortex equations is that the dynamical timescales of the 
disturbances are long compared to the local disk rotation 
time. This allows for an analysis which can focus on the long 
term properties of vortices under these circumstances free 
of the effect of soundwaves. Soundwaves and other short 



Table 1. Comparison of the various scalings employed that 
result in some of the approximations used in the literature 
to study small sections of accretion disks. 



Quantity 


SSB 


LSB 


Elongated 
Vortices 


Anelastic 
Model 


This 
Work 


e = c s /v k 


< 1 


« 1 


« 1 


<C 1 


< 1 


o . Maen # 


<C £ 


-1 

~ 1 


e <c 8 <c 1 


e <C 8 <C 1 


~ 1 


azimuthal 
scales 


oHa) 


r T> 

eKq 


eKo 


oeHo 


r> 

-no 


vertical 
scales 


z r> 

oHo 


7") 

eHo 


eito 


£ Ti 

oeHo 


eKo 


radial 
scales 


otto 


eito 


oerio 


Ac F?~ 


eito 


azimuthal 


fir 
< ><- s 


Cs 


fir 


fir 




vertical 
speeds 


Sc s 


c n 


5c s 


8c s 


fC - 


radial 
speeds 


Sc s 


c s 


5 2 c s 


Sc s 




temporal 
scales 


i/n 




i/(sn ) 


i/n 


i/(en ) 


dynamical 
pressure 


Spcl 


pel 


5 2 pcl 


8 2 pcl 


pel 


dynamical 
density 


P 


p 


8 2 P 


S 2 P 


P 



time fluctuations in this study have been filtered out by the 
fact that we are considering time scales which are longer 
than what typifies these wave phenomena. It is this long 
timescale view together with the assumption made here 
that the radial and vertical velocity fluctuations are smaller 
than the azimuthal ones that ultimately lead to the shallow 
water theory limit arrived at here. 

6.2. Walls, defects, edgewaves and instability 

In an examination of single PV defect disturbances of an 
infinite domain we saw that by letting the PV on one side of 
the defect tend to infinity is mathematically equivalent to 
replacing the defect with an impenetrable wall. Indeed, the 
limiting form of the Rossby edgewave wavespeed (|47[) in the 
limit where the PV on one side becomes infinite recovers the 
expression for the wavespeed on the semi-infinite domain 
([55]) . It would seem then, at least within the context of this 
theory, that this equivalency may be a useful conceptual 
tool. For instance, the use of impenetrable boundary condi- 
tions is widely considered to be artificial in most cases other 
than for terrestrial laboratory flows in a channel. This con- 
sensus is especially true in the study of astrophysical disks. 
Of course as disks no not have delineated solid boundaries, 
this view is physically justified. However, it is conceivable 
that disks may have regions in which the radially varying 
potential vorticity undergoes rapid changes. One obvious 
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candidate region would be the star-disk boundary. If the 
star rotates significantly sub-Keplerian then the PV in the 
boundary layer region is likely to be significantly different 
than the PV immediately outside the boundary layer. The 
transition between the two can then be viewed as a de- 
fect in the disk's PV and in the limit where the PV in the 
boundary layer greatly exceeds the PV outside, the argu- 
ments laid out here could be used to justify replacing the 
defect with a solid boundary for the purposes of analyzing 
the response in the extended disk. 

With this in mind, we can view the edgewave prop- 
agating along the solid boundary as a limiting form of 
a Rossby edgewave propagating along a defect. It follows 
that much of the machinery involved with the principle of 
counterpropagating Rossby waves (CRW), originally pro- 
posed in Hoskins et al. (1985) in the context of understand- 
ing barotropic/baroclinic instabilities in planetary atmo- 
spheres, should carry over to understanding disks as well. 
In Section 3 we presented a normal mode analysis on a finite 
domain bounded by two impenetrable walls. The instabil- 
ity associated with it is a manifestation of the SRI. In view 
of the analysis, we see that the two modes appearing un- 
der those circumstances are Rossby edgewaves propagating 
along the boundary walls. Under the right conditions (i.e. 
wall separation) these edgewaves interact with each other in 
such a way as to promote an instability in the way outlined 
and demonstrated in Molcmaker et al. (2001) - much in the 
same way that edgewaves are responsible for the emergence 
of the PPI (Goldreich et al., 1985, and see below). 

Now, given the foregoing argument, the edgewaves as- 
sociated with the boundaries can be viewed as a limiting 
form of Rossby edgewaves propagating along a PV defect 
in which the PV on one side greatly exceeds the other. 
In this way, the left boundary looks like it has a negative 
PV gradient while the gradient on the right boundary looks 
positive which satisfies the two-signed PV-gradient instabil- 
ity condition of Charncy & Stern (1962) which means that 
the mechanism of the instability can be physically rational- 
ized as arising from the interaction of CRWs (Hoskins et 
al., 1985, Heifetz et al., 2004). Such instabilities ultimately 
draw their energy from the background shear - without 
which the instability vanishes as is the case here. An ex- 
amination of the SRI in which the channel boundaries are 
replaced with either PV defects or strongly varying smooth 
PV profiles is the subject of a future study. 

A real disk probably is rotationally supported all the 
way through, which would mean that the only way that 
one could realistically argue for significant changes in its 
PV profile would be if the steady-state height of the disk 
changes as a function of the radial coordinate. Such local- 
ized annuli will probably introduce some departures from 
Keplerian flow as well. The possibility of this scenario has 
been proposed before by Lovelace et al. (1999) and Li et 
al. (2000) in the context of two-dimensional disks charac- 
terized by strong radial entropy gradients where they show 
that Rossby waves turn unstable under suitable circum- 
stances. It is in the opinion of this author that the viability 
of this process shows promise and ought to be further con- 
sidered. 

6.3. Dynamical results in comparison to other studies 

The linear instability discovered by Papaloizou & Pringlc 
(1984,1985) was analytically demonstrated in an analysis of 



homcntropic gas in slender accretion tori (i.e. not exactly 
rotationally supported) constrained by pressure conditions 
on the boundaries of the domain. This means to say that 
the vertical and radial scales are much smaller than the 
azimuthal scales and the flow is contained in a narrow an- 
nular channel. Furthermore, the modal disturbances have 
very nearly uniform potential vorticity throughout and have 
the character of being essentially vertically hydrostatic as 
illucidated by Goldreich et al (1986). This instability oc- 
curs for sinuous modes (i.e. having even symmetry with 
respect to the disk midplane) while no such instability is 
observed for varicose modes (i.e. having odd with respect 
to the disk midplane). The instability has been physically 
rationalized by understanding that the pressure boundary 
conditions applied at the inner and outer walls of the undu- 
lating tori essentially set up a situation in which edgewaves 
countcrpropagate with respect to each other and interact 
across a region which is evanescent - leading ultimately to 
instability (Goldreich et al, 1986, Blaes & Glatzel, 1986). In 
view of this mechanistic understanding, the PPI should be 
considered the disk analog of Eady edgewaves encountered 
frequently in terrestrial atmospheric studies (e.g. Davies & 
Bishop, 1994). 

The dynamics considered here are for rotationally sup- 
ported disks (q = 3/2) which differ mathematically from 
tori (q ^ 3/2) which have significant radial pressure support 
at the lowest order of the dynamics. The instability of uni- 
form potential vorticity modes discussed in Section 3 and 
which is an incarnation of the strato-rotational instability 
(SRI, Dubrullc, 2004) is not the same as the Papaloizou- 
Pringle instability (PPI). The main differences are (i) the 
symmetry of the modes discussed in this work are varicose 
while they are sinuous for the PPI and (ii) the PPI arises 
in part from the pressure boundary conditions at the inner 
and outer boundaries while the instability here arises from 
the no normal flow condition at the inner and outer bound- 
aries. Both instabilities, however, share (i) the character- 
istic as arising from the interaction of counterpropagating 
edgewaves of some sort and are explained in terms of the 
Hayashi- Young criterion for wave instability in shear flow 
(Hayashi & Young, 1987) and, (ii) the character of belong- 
ing to a general class of barotropic/baroclinic instabilities 
related to the dynamics of potential vorticity disturbances 
(Brctherton, 1966, Hoskins et al., 1985, Davies & Bishop, 
1994). 

6.4. Future directions 

A natural test and evaluation of the shallow water the- 
ory developed in this study will be to study the nonlinear 
solutions of non-uniform potential vorticity disturbances. 
At the outset such an investigation will surely require nu- 
merical computational simulations. The question of what 
boundary conditions to use will be a sensitive one since, for 
instance, certain conditions may lead to strong linear in- 
stabilities and the relevance of them to real disks are ques- 
tionable (although certainly not yet settled). For example, 
by requiring no-normal flow boundary conditions can excite 
disturbances which include uniform potential vorticity per- 
turbations which have been shown to be the analog of the 
SRI. However, in light of the arguments presented here, it is 
possible that such boundary conditions and their resulting 
dynamical response could be qualitatively representative of 
what goes on in bonafide real disk systems. 
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It should also be considered that the shallow water the- 
ory developed here (by its construction) is unable to follow 
the dynamical response of disturbances in which there is 
a vertical variation in the entropy on account of the as- 
sumptions made about the thermodynamical response of 
the fluid made at the outset (i.e. the assumption of a poly- 
tropic gas). Vertical variations of this sort may be playing 
a critical role in the development of the dynamics reported 
in Barranco & Marcus (2005) since it was observed in that 
work that vertically propagating disturbances contribute 
to the complex vortex dynamics seen to manifest in their 
simulations ( see also the discussion on the stability of tall 
columnar vortices in Dritschel & Juarez, 1996). As one of 
the long term goals we have in mind for developing a shal- 
low water formalism of non-axisymmetric disturbances is 
to better understand the emergence of such interesting dy- 
namics, it will be worthwhile to extend the theory devel- 
oped here in order to accommodate such circumstances by 
allowing for vertical variations of some type beyond the fluc- 
tuations of the disk height we follow here. One way would 
be to develop a multi-stacked shallow water theory in a 
manner similar to what has been developed in the study of 
MHD disturbances in the solar tachoclinc (Gilman, 2000, 
Cho et al, 2007). 

We would like to make one final observation pointed out 
by the referee of this article. With respect to the problem 
of stellar tachoclines Zaqarashvili et al. (2007) discuss how 
the inclusion of a background azimuthal magnetic field can 
split an ordinary Rossby wave into two different parallel 
modes. By analogy a similar phenomenon is likely to occur 
in magnetized disks (see also Matsumoto & Tajima, 1995). 
A strategy like the one employed in this work may be used 
to further detail the scope and shape of the modification of 
Rossby waves in such circumstances. 
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